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Abstract 

We study effective energies of heavy-light meson correlation functions in lattice QCD 
and a small volume of (0.2 fm) 4 to non-perturbatively calculate their dependence on the 
heavy quark mass in the continuum limit. Our quenched results obtained here constitute 
an essential intermediate step of a first fully non-perturbative computation of the b-quark's 
mass in the static approximation that has recently been presented as an application of a new 
proposal to non-perturbatively renormalize the Heavy Quark Effective Theory. The renor- 
malization constant and the improvement coefficients relating the renormalized current and 
subtracted quark mass are determined in the relevant parameter region at weak couplings, 
which allows to perform the numerical simulations at several, precisely fixed values of the 
renormalization group invariant heavy quark mass in a range from 3 GeV to 15 GeV. 
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1 Introduction 



During the last years, the elementary particle physics community has seen a growing interest 
and activity in the study of B-meson physics and its phenomenology. While on the experi- 
mental side this interest reflects in the various facilities that are currently running to explore 
CP-violation in the B-system [1-3], it is nourished on the theoretical side by the demand to 
determine transition matrix elements of the effective weak Hamiltonian in order to interpret 
the experimental observations within (or beyond) the standard model and to (over-)constrain 
the unitarity triangle. For the computations of such matrix elements between low-energy 
hadron states to become valuable contributions in this field, they have to be carried out 
non-perturbatively, which is the domain of QCD on the lattice. However, in contrast to light 
quarks which as widely spread objects are predominantly exposed to large-volume limita- 
tions, heavy quarks are extremely localized ~ 1/(4 GeV) — 0.04fm) and thus also 
require very fine lattice resolutions, because otherwise one would face huge discretization er- 
rors. Mainly for this reason, realistic simulations of heavy-light systems involving a b-quark 
(even in the quenched approximation) are impossible so far [4]. 

A theoretically very appealing way out of this restriction is to recourse to the Heavy 
Quark Effective Theory (HQET) [5,6]. This comes at a prize, though. As a consequence of 
its different renormalization properties, physical quantities deriving from expectation values 
calculated in the effective theory are affected by power-law divergences in the lattice spacing 
that can not be subtracted perturbatively in a clean way: the continuum limit ceases to 
exist, unless the theory is renormalized non-perturbatively [7]. 

Only recently, an approach to overcome these deficiencies has been developed from 
a solution for the problem of a completely non-perturbative computation of the b-quark 
mass in the static approximation including the power-divergent subtraction [8] to a new 
method addressing the general class of renormalization problems in HQET [9]. At its root 
stands the insight that these power-law divergences can be removed by a non-perturbative 
matching procedure to relativistic QCD in a finite volume. In fact, the smallness of the 
physical volume represents a characteristic feature in this strategy — since only then one is 
capable to incorporate and simulate the b-quark as a relativistic fermion — and motivates 
the investigation of QCD observables in a small-volume setup presented in this work. Here 
we will concentrate on the non-perturbative heavy quark mass dependence of effective heavy- 
light meson energies in the continuum limit, whose numerical knowledge is crucial to apply 
the proposal of Ref. [9] to the determination of the b-quark's mass in leading order of HQET. 1 
An extension to also examine the mass dependence of a few more quantities, which aims at 
quantitative non-perturbative tests of HQET by comparing static results with those obtained 
along the large quark mass limit in small- volume QCD, is in progress [12]. 

In Section [21 we first recall the main ideas of the general matching strategy between 
HQET and QCD of Ref. [9] (Section 12. 1|) . then introduce our observables considered (Sec- 
tion and finally describe how these can be calculated as functions of the (renormalized) 
mass of the heavy quark by numerical simulations (Section 12.3(1 . Section 13.11 is devoted to 
some intermediate results on a renormalization constant and improvement coefficients that 
are needed to renormalize the heavy quark mass in the relevant parameter range. Our central 
results on the mass dependence of the heavy- light meson energies are discussed in Section T3. 21 

1 For another method to determine and Fb, which also starts from lattice QCD in small volume but 
employs extrapolations of finite- volume effects in the heavy quark mass, see Refs. [10,11]. 
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and we conclude in Section 0] 

In the calculations reported here we still stay in the quenched approximation. We want 
to emphasize, however, that this study of the heavy quark mass dependence of suitable 
observables (as an important part of the general non-perturbative approach to HQET of [9]) 
can be expected to be numerically implementable as well for QCD with dynamical fermions 
at a tolerable computing expense, because on basis of the experiences made in Refs. [13, 14] 
the use of the QCD Schrodinger functional brings us in a favourable position where in 
physically small to intermediate volumes dynamical simulations are definitely easier than 
with the standard formulation. 

2 Computational strategy 

For the rest of this paper we will focus on meson observables derived from heavy-light 
correlation functions in finite volume and their dependence on the heavy quark mass in the 
continuum limit. Before we explain in detail how this dependence can be singled out in an 
actual numerical computation, we want to clarify the special role of QCD in a finite volume as 
a material component of the more general idea advocated in Ref. [9] to non-perturbatively 
renormalize HQET; therefore, it suggests itself to briefly summarize this idea in the first 
subsection. 

2.1 From QCD in small volume to non-perturbatively renormalized HQET 

A long-standing problem with lattice computations in HQET is the occurrence of power-law 
divergences during the renormalization process (cf. Refs. [15-17]), implied by the allowed 
mixings of operators of different dimensions coming with coefficients {c^} that contain inverse 
powers of the lattice spacing a. At each order of the HQET expansion parameter (1/m, 
where m is the heavy quark mass), new such free parameters Cfc arise, which in principle are 
adjustable by a matching to QCD; but owing to incomplete cancellations when performing 
this matching only in perturbation theory, one is always left with a perturbative remainder 
that still stays divergent as a — > 0. Therefore, the continuum limit does not exist. 

Already HQET in leading order, the static approximation, exhibits this unwanted phe- 
nomenon. In this case the kinetic and the mass terms in the static action mix under renor- 
malization and give rise to a local mass counterterm 8m ex 1/a, the self-energy of the static 
quark, which causes a linearly divergent truncation error if one relies on an only perturbative 
subtraction of this divergence. A prominent example for a quantity suffering from it is the 
b-quark mass itself: past computations in the static approximation [18] were limited to finite 
lattice spacings, and the continuum limit was impossible to reach. 

A viable strategy to solve this severe problem of power divergences is provided by a non- 
perturbative matching of HQET and QCD in finite volume [8,9]. To integrate the present 
work into the broader context of Ref. [9], we reproduce the central line of reasoning here. Let 
us consider QCD consisting of (generically iVf — 1) light quarks and a heavy quark, typically 
the b-quark. In the effective theory, the dynamics of the heavy quark is governed by the 
HQET action, which formally is an expansion in inverse powers of the heavy quark mass. 
(For further details see e.g. Section 2 of Ref. [9].) On the renormalized level, the effective 
theory discretized on a lattice can be defined in terms of parameters {ct} that comprise 
those specifying the light quark sector (e.g. the bare gauge coupling, g$ = ci, and the masses 
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of the light quarks) and coefficients of local composite fields in the 1/m-expansion of the 
lattice action, supplemented by further coefficients belonging to local composite operators 
which will be needed when including their correlation functions into the renormalization 
program. In other words, if these parameters are chosen correctly, HQET and QCD are 
expected to be equivalent in the sense that $ H Q ET (M) = $Q CD (M) + 0(1/M n+1 ) holds for 
suitable observables $ in both theories, where for simplicity only the dependence on the heavy 
quark mass, here represented by the scheme and scale independent (and thus theoretically 
unambiguous) renormalization group invariant quark mass, M, is stressed. Now suppose 
that in some way the parameters of QCD have already been fixed to proper values. Then 
the parameters {c^} in the effective theory may just be determined through its relation to 
QCD by requiring a set of matching conditions: 

<D* QET (L,M) = ^ CD (L,M), k = l,...,N n . (2.1) 

In this equation, the index k labels the elements of the parameter set {c^} defining the 
effective theory up to l/M n+1 -corrections (where, for instance, the additive mass renormal- 
ization 8m mentioned above is amongst them), and the conditions (|2.1j) determine the Ck 
for any value of the lattice spacing. Moreover, we have already indicated the dependence 
of the observables <3?fc in eq. (|2.1jl on another variable which will become important in the 
following: the linear extent L of a finite volume. 

To substantiate this L-dependence, we note that in order to circumvent the difficulties 
with the power-law divergences from the start, the matching equations (|2.1|) are understood 
as non-perturbative conditions in which both sides are to be calculated with the aid of 
numerical simulations. From the practical point of view, this in turn also means that one 
must be able to simulate the b-quark as a relativistic fermion. Hence, the linear extent L of 
the matching volume (i.e. where eq. (|2.1|) should hold) has to be chosen with care. On the 
one hand, it should fulfill LM 3> 1 to apply HQET quantitatively on the l.h.s. of (|2.1jl . i.e. to 
suppress 1/m-corrections, and on the other hand one has to ensure aM < 1 on the r.h.s. to 
treat the heavy quark flavour in the relativistic theory and avoid large lattice artifacts so 
that the continuum limit is under control. Taking these constraints together, while at the 
same time keeping the number of lattice points manageable for present-day computers, one 
then ends up with a volume for imposing the matching conditions (|2.1|) that is physically 
small; in our application later it will turn out to be of the order of (0.2 fm) 4 . 

Having highlighted QCD in finite volume as key ingredient for the practical realization of 
the non-perturbative matching, we close this subsection with a few remarks on the subsequent 
(but not less important) steps that are involved to adopt it as a general approach for non- 
perturbative computations using the lattice regularized HQET. These steps, together with 
two applications of the proposal as a whole, are also worked out in Ref. [9], which the reader 
should consult for a thorough discussion. 

• Rather than directly identifying the quantities entering Q2.1j) with physical, experi- 
mentally accessible observables (by which one would sacrifice the predictability of the 
effective theory), they must be properly chosen as renormalized quantities computable 
in the continuum limit of lattice QCD, which in turn necessitates the use of a small vol- 
ume as explained before. Of course, their explicit form still depends on the application 
in question. One may think e.g. of hadronic matrix elements or, as in the computa- 
tions reported in the following sections, effective masses (respectively, energies) that 
are deduced from the decay of two-point heavy- light correlation functions. 
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• Apparently, although these matching conditions to fix the parameters of HQET are 
perfectly legitimate (since the underlying Lagrangian does not 'know' anything about 
the finite volume), one still has to make contact to a physical situation, where the inter- 
esting quantities of the infinite- volume theory can be extracted at the end. Employing 
the same lattice resolution as in the small volume computation, however, would again 
demand too many lattice points. This gap between the small volume with its fine 
lattice resolution on the one side and larger lattice spacings, and thereby also larger 
physical volumes, on the other is bridged by a recursive finite-size scaling procedure 
inspired by Ref. [19]. As a result, the & k are obtained at some larger volume of ex- 
tent L = 0(1 fm), where the resolutions a/L are such that at the same lattice spacings 
(i.e. at the same bare parameters) volumes with L ~ 2fm — to accommodate physical 
observables in the infinite-volume theory — are affordable. 

• Finally, the approach requires a physical, dimensionful input. This usually amounts 
to relate the observables <3? fc of the effective theory calculated in large volume to 
some experimental quantity. Which quantity this actually might be has to be decided 
when a concrete application is addressed. (E.g., in the application to compute %, it 
is the B-meson mass [9].) 



2.2 Lattice setup and observables 

In our investigation of QCD in a small volume we distinguish between a light (T) and a 
heavy ('h') quark flavour. The lattice regularized theory is formulated in a Schrodinger 
functional (SF) cylinder of extent L 3 xT [20,21]: the gluon and quark fields are subject to 
periodic (Dirichlet) boundary conditions in spatial (temporal) directions, and — if not ex- 
plicitly stated otherwise — we assume T = L from now on. In principle, the aforementioned 
matching strategy between HQET and QCD is not restricted to the SF as the only possi- 
ble finite-volume scheme to treat the involved heavy-light systems in the relativistic theory, 
but in this way its practical implementation will much profit from known non-perturbative 
results on the renormalized quantities that will be needed in the following. 
Starting from the 0(a) improved heavy- light axial and vector currents, 

(Aij^x) = ^(x^jsip^x) + ac A ^(d,j, + dl) {ip } {x)j 5 i/) h (x)} , (2.2) 
(Vi)fj,(x) = ip^x^^tp^x) +acy~(d u + 8%) {iip^a^ iph(x)} , (2.3) 

(where numerical values for the coefficients ca and cy in the quenched approximation are 
taken from Refs. [22,23]), we construct their correlation functions in the SF [24,25] as 

fA(xo) = -yE<(^)°(*)Ch(y)75Ci(z)>, (2.4) 

kv(xo) = -^^<(Vi) fc (x)C h (y) 7fe Ci(z)>, (2.5) 

schematically drawn in Figure ^ Based on these, we define L-dependent energies in both 
channels by combining the usual forward and backward difference operators to take the 
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x = 



x = L 



Figure 1: Illustration of the correlation function /a, while in ky the insertion of (-Ai)o is replaced by (Vi)k- 

logarithmic derivatives 

r PS (L,M) = -l(do + d* )ln[f A (x ,M)} , (2.6) 

A xo=L/2 

r v (L,M) = -hdo + d* )ln[k Y (x ,M)} T = L, (2.7) 

A xq=L/2 

in which all multiplicative renormalization factors drop out. These observables also de- 
pend on the bare coupling, go, and — as the light quark mass is put to zero in the actual 
computations — on the heavy quark mass variable which, as already in Section 12.11 one 
may conveniently set to be the renormalization group invariant (RGI) mass of the heavy 
quark, M. Having in mind the corresponding correlators in the effective theory, where in 
the static (i.e. infinite-mass) limit the two currents fall together owing to the heavy quark 
spin-symmetry, we also form their spin-averaged combination 

r av = \ (r PS + 3r v ) , (2.8) 

which is expected to deviate from the static limit by the smallest 1/m-corrections [9]. 



2.3 Determination of the heavy quark mass dependence 

As mentioned before and detailed in Refs. [8,9], the non-perturbative matching between 
HQET and QCD can be applied to determine the mass of the b-quark in the static approxi- 
mation. To achieve this goal, the quantities in eq. (|2.1|) have to be identified with suitable 
observables that can be calculated by numerical simulations. Since we want to evaluate them 
in the continuum limit, we have to ensure a situation with fixed renormalized parameters, 
along which the limit a/L — > may be performed. In regard of the (suppressed as well as 
indicated) dependences in the matching conditions (|2.1|) . this means that for each value of 
L/a we have to make particular choices for the gauge coupling and the quark masses in the 
light and heavy sectors, respectively. The first two of them are immediately offered by 

^hqet = - 2 ( L ) = congtanti (2.9) 
^HQET = m[ = Q) ( 21Q ) 

where g 2 denotes the renormalized finite-volume coupling in the SF scheme that runs with 
the box size L [26] and m\ = mf CAC is the PCAC current quark mass of the light flavour. 
These conditions allow to readily benefit from earlier work, because when fixing g 2 to one of 
the values used for the (quenched) computation of the non-perturbative quark mass renor- 



malization in [27], they translate into triples (L/a,f3 = 6/<? 2 ,/-ti) that can directly be taken 
over from this reference. 2 

For the computation of the b-quark's mass in leading order of HQET one yet needs one 
more condition in order to fix the parameter a5m in the static Lagrangian. Following [8,9], 
the most natural candidates for this purpose are now the pseudoscalar and spin-averaged 
energies of the form (|2.6|) and (|2.8[) introduced in Section T2. 21 For k = 3, eq. (|2,1|) then turns 
into 

L (r stat + m) = f ™(L, M) = <^ CD (L, M) = L T x , (2.11) 

where Tx = Tx(L, M), X = PS, av, are defined in QCD and small volume with a relativistic 
b-quark and carry the entire quark mass dependence. On the l.h.s., r sta t denotes the analogue 
of Tx in the effective theory, where the heavy quark entering the correlator (|2.4j) is treated 
in the static approximation, and implicitly contains the problematic, linearly divergent mass 
counterterm abm. The matching equation (|2.11j) may be connected to a physical quantity in 
large volume (the B-meson mass [9]) by finite-size scaling (cf. Section E^). But without going 
into further, unnecessary details we only note that this yields a (dimensionless) equation, 
which apart from Tx(L, M) only contains energy differences computable in the static theory 
— such that the counterterm a5m cancels out — and the experimental mass of the B-meson. 
Since all pieces in that equation can be evaluated non-perturbatively and in the continuum 
limit, the quantitative knowledge of Tx(L, M) as a function of M in the relevant quark mass 
region finally allows to solve it for the desired numerical value of the non-perturbatively 
renormalized RGI mass of the b-quark, Mb, in the static approximation [9]. 

So the task is here to non-perturbatively determine the heavy quark mass dependence 
of the effective heavy-light meson energies Tx{L,M). It thus remains to discuss how we fix 
the linear size of the finite volume, L = Lq, and a set of dimensionless quark mass values, 

z = L M = constant , (2.12) 

at which the numerical computation of Tx for various lattice spacings takes place in practice. 

We answer this question again on basis of the present numerical knowledge in the frame- 
work of the SF. Among the available constant values of the renormalized SF coupling that 
belong to different sets of triples (L/a, f3, k\) known from Ref. [27] we fix, as was initially 
required by the condition 1)2. 9|) . 

f(L /2) = 1.8811, (2.13) 

where the (inverse) lattice spacing varies within 6 < Lq/(2o) < 16. This choice is particularly 
convenient, because from the same reference one then finds that Lq is given in units of a 
certain maximal size, L max , of the SF box, namely 

L = = o.36 r « 0.18 f m , (2.14) 

which is implicitly defined by g 2 (L max ) = 3.48 and, through its relation to the hadronic 
radius tq = 0.5 fm [28], also known in physical units: L max /ro = 0.718(16) [29]. Indeed, from 
a closer inspection of the parameters corresponding to the condition (|2.13|) (see Table ^) 

2 In this case the lattice action's hopping parameter associated with the light quark flavour, k\, just equals 
the critical hopping parameter, k c , as a consequence its very definition: mi(/c c ) = 0, see Ref. [27]. 



6 



one may convince oneself that an accompanying volume of some small extent as (|2.14|) is a 
sensible compromise between the associated resolutions being small enough to comfortably 
accommodate a propagating (albeit the system size confining) b-quark and being able to 
consider heavy quark masses large enough to suppress 0(l/z)-corrections (i.e. z ^> 1 but 
not too large to avoid a substantial enhancement of cutoff effects) as well as to cover the 
physical RGI b-quark mass scale itself. 

In order to satisfy eq. (|2.12|) for the intended numerical values z while L/a and (3 are 
changed to approach the continuum limit, we still have to relate the RGI mass of the heavy 
quark, M, to the corresponding hopping parameter as the simulation parameter for the 
bare heavy quark mass in the lattice action. One possibility would be to employ the (0(a) 
improved) relation between the RGI and the bare PCAC quark mass, m^, viz. 

\* M Z A (g )(l + b A am qM ) , 2 . 

M = — u — rxm h + 0(a), fi = 1 L , (2.15) 

m h (/i ) Z P (g , L )(l + bpam^ h ) 

where mh(^o) denotes the running heavy quark mass renormalized at the scale /io = l/-^o 
in the SF scheme and (the inverse of) its ratio to the RGI mass M, the first factor on 
the r.h.s., is inferable from the literature as will be elaborated on after eq. (|2.19|) . The 
numerator and denominator in the second factor account for the renormalization of the axial 
current and the pseudoscalar density in the 0(a) improved theory, respectively. Here, the 
subtracted bare heavy quark mass is defined as usual by m q ^ = tuq ^ — m c with m c the 
critical value of the bare quark mass mo,h, which in terms of the hopping parameters is given 
by am qj h = — k^ 1 ). Although the involved renormalization factors are in principle 

known from Refs. [27,30] in the quenched case, this would demand additional simulations to 
appropriately tune — and thereby the PCAC mass as a secondary quantity composed 
of correlation functions (cf. eq. Q3.3|) ) — until the condition ()2.12j) is met with sufficient 
precision. 

However, a safer (and also more economic) way to estimate such that a constant 
z = LqM can be enforced in advance without any additional tuning runs is as follows. Recall 
that alternatively the renormalized quark mass can be written in terms of the 0(a) improved 
subtracted bare quark mass, m qj h, as 

fn h = Z m m qjh , m qjh = m qjh (l + b m am^ h ) , (2.16) 

where non-perturbative estimates on the improvement coefficient 6 m and the finite combina- 
tion of renormalization constants 

Za(#o) 

have been published in Ref. [30] for the quenched approximation. The latter identity may 
then be used to eliminate the (in the SF scheme unknown) renormalization factor Z m in 
eq. (|2.16jl in favour of numerically known ones. The decomposition of the RGI heavy quark 
mass analogous to ()2.15(l now assumes the form: 

M = Z( f^ Z ^ x m q , h (l + fe m am q , h ) + O (a 2 ) , Mo = 1/L • (2-18) 

m h (/i ) Zp(g ,L ) 

The last piece to be addressed in this equation is the universal, regularization independent 
ratio of the RGI quark mass to the renormalized mass at fixed renormalization scale, which 
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for later reference we call 



M 

h(L ) = — — , mo = l/L . (2.19) 

In a mass independent renormalization scheme such as the SF it is flavour independent and, 
according to Ref. [27], can be expressed as 



z P {L ) r _ 2 ^-^0/(260) / r w , 



(2.20) 



with Lq = L max /2, /i = 2 6 /L max and &o>^o the leading-order coefficients in the perturbative 
expansions of the renormalization group functions of the running coupling and quark mass, 
(3{g) and r{g). In the ratio of Zp-factors we encounter a scale evolution through changes by 
finite step sizes in the SF renormalization scale L = l/ju that has been non-perturbatively 
computed in [27], whereas the remainder lives in the high-energy regime where the coupling 
is reasonably small to evaluate it in perturbation theory. Since in eq. (|2.14jl we just chose the 
linear extent Lq of the small volume, where the non-perturbative matching to the relativistic 
theory is to be performed, as a proper multiple of the scale L max , a numerical estimate for 
h(Lo) in the continuum limit, eq. (|2,19|) . can be directly extracted from the results already 
published in [27] and will be quoted in Section f3.2l later on. 

Our simulation parameters, which fix via eqs. (|2.9|) . (12.10(1 and (|2.13j) the physics in 
the relativistic sector and via eq. (|2.14|) the size of the matching volume where to determine 
the non-perturbative heavy quark mass dependence of our observables, are summarized in 
Tabled Values for the other quantities, which have to be inserted into eq. (|2.18|) in order to 



set 


L/a 


P = 6/sg 




t(L/2) 






A 


12 


7.4082 


0.133961(8) 


1.8811(22) 


2.397(17) 


0.6764(6) 


B 


16 


7.6547 


0.133632(6) 


1.8811(28) 


2.393(18) 


0.6713(8) 


C 


20 


7.8439 


0.133373(2) 


1.8811(22) 


2.379(22) 


0.6679(8) 


D 


24 


7.9993 


0.133159(4) 


1.8811(38) 


2.411(20) 


0.6632(8) 


E 


32 


8.2415 


0.132847(3) 


1.8811(99) 


2.397(52) 


0.6575(13) 



Table 1: Our parameter sets that refer to the light quark sector and have fixed SF coupling, g 2 (L/2) = 1.8811. 
The parameters L/a, f3 and k\ = k c of A, B, D and E are taken over from [27] without changes, whereas 
we performed new simulations to add with C a further lattice resolution for this work. Moreover, the 
renormalization constants Zp — Zp(go, L/a) differ from those of [27] in that they have been computed 
in the context of [31] using the two-loop value [32] for the boundary improvement coefficient ct. 

obey the condition (|2.12|) of constant z = LqM, will be specified when we come to describe 
the actual calculation of the z-dependence of Fx = Tx{Lo, M), X = PS, av, in Section I5~2l 
At this stage we still want to direct the reader's attention to an issue that can already 
be foreseen to become a potential limitation in any application of the numerically computed 
quark mass dependence: even with eq. (|2,18|) . the dimensionless mass parameter z can only 
be fixed through quantities known to a certain precision. In particular, the improvement 
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coefficient b m and the renormalization factor Z are only available from Ref. [30] 3 for 6.0 < 
(3 < 6.756, which are significantly lower than the /3-s of Table ^ we are interested in. It 
hence appears quite difficult to reliably guess from the numbers in [30] values for b m and Z 
in the /3-region relevant here and, even more, to quantitatively assess the additional error 
contribution those estimates would then be afflicted with. Since this constitutes a dominant 
source of uncertainty, which would also propagate into any quantity that explicitly derives 
from the z-dependence of the observables studied here — such as, for instance, the b-quark 
mass through the matching of HQET and QCD [9] — , we first determine 6 m and Z exactly 
for the /3-values of Table Q and also improve their numerical precision. 

3 Results 

We now present the numerical results of our (quenched) computation of the heavy-light 
meson energies Tx(Lq, M), X = PS,av, which — as outlined in the previous subsection - 
basically consists of two parts: first, we determine b m and Z in the /3-range relevant for a 
matching in physically small volume in order to be able to fix the RGI heavy quark mass 
to a set of desired values z = LqM with sufficient precision, and second, the LqTx(Lo, M) 
are calculated at these values of z in the continuum limit. This will then allow for smooth 
representations of the energies as functions of z, which eventually may be used to interpolate 
them to the b-quark scale, z^ = LqM^ [9]. 

3.1 Coefficients 6a — bp, b m and renormalization factor Z for 7.4 < (3 < 8.2 

Here we proceed in the spirit of Ref. [30] , where the idea of imposing improvement conditions 
at constant physics was advocated. In that work, which considered a range 6.0 < [3 < 6.756 
of bare couplings commonly used in large- volume simulations, this was realized by keeping 
constant the ratios L/r$ and T/L = 3/2 (supplemented by the SF-specific choices C = C = 
and 9 = 0.5). However, despite ln(a/ro) as a function of (3 needed to fix L/tq to some 
suitable value for given L/a is available for 5.7 < (3 < 6.92 [35] and even beyond [36], the 
condition of constant L/r$ can not be transfered to the present situation. The reason lies 
in the fact that for /3-values in the range we are interested in, 7.4 < (3 < 8.2, enforcing 
an improvement condition such as L/tq = 0(1) would lead to prohibitively large values 
of L/a in the simulations. Therefore, to replace the latter, we exploit one of the already 
built in elements of the matching strategy between HQET and QCD as explained in the 
foregoing section: namely, as for the computation of the energies Tx on the QCD side we 
have to work along a line of constant physics in bare parameter space anyway, we can directly 
adopt the pairs (L/a, (3) of Tabled which correspond to L/L max = 1/2 = 0.36ro/L m ax and 
simultaneously to a constant renormalized coupling of g 2 (L/2) = 1.8811 (see eqs. (|2.14|) and 
(|2.13[0 . Within the present application, this constitutes a much more natural and equally 
admissible choice of improvement condition and thereby, in the same way as in Ref. [30] , the 
improvement coefficients 6a — bp and b m as well as the renormalization constant Z become 
smooth functions of g$ but exactly in the region where they are needed. 

3 The /3-range considered in [30] was chosen to cover values commonly used to simulate O(a) improved 
quenched QCD in physically large volumes. Similarly, the results from determinations through chiral Ward 
identities with mass non-degenerate quarks [33,34] refer to this region and are not applicable here, too. 
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Following [30] — and also referring to this reference for any unexplained details — , the 
improvement coefficient 6 m and the renormalization constant Z (as well as the difference of 
coefficients 6 a — bp, though it does not enter the subsequent computations in Section 13 .2[) 
can be determined by studying QCD with non-degenerate quarks. This is particularly ad- 
vantageous in case of the quenched approximation at hand, since then the structure of the 
0(a) improved theory stays quite simple. For instance, the improvement of the off-diagonal 
bilinear fields = X 1 ± i X 2 , X = A^, P, emerging as a consequence of the broken isospin 
symmetry, is the same as in the degenerate case, except that the 6-coefficients now multi- 
ply cutoff effects proportional to the average ^(am qi j + amy) of the subtracted bare quark 
masses, m qi j = mo,j — m c , which themselves are separately improved for each quark flavour: 

m qii = m qii (l -(- b m am qi i) . (3.1) 

(Here and below the indices i,j label the different quark flavours.) Identifying, for instance, 
the flavours in the isospin doublet with a light and a heavy quark as before, the corresponding 
PCAC relation reads 

d^A±(x) = (mi + m h )P ± (x) , (3.2) 

and the renormalization constants Z\ and Zp that come into play upon renormalization are 
just those known in the theory with two mass degenerate quarks. 

Accordingly, the fermionic correlation functions defined in the SF and involving the 
axial current and the pseudoscalar density [24] generalize to f^(xo) = — \ (Aq(x)0~) and 
fp (xo) = — \ (P + (x)0~), with pseudoscalar boundary sources decomposed as = O 1 ± 
iO 2 where O a = o 6 S yz C(y)75 ^ ra C( z )- Then the improved bare PCAC (current) quark 
masses 4 as functions of the timeslice location xq are given by 



, .,. fr ... d f/{x ) + ac A dodofp J (x ) 
mijixo-L/a, T/L, 9) = , (3.3) 

2/p / (a;o) 



where only here we explicitly indicate their additional dependence on L/a, T/L and the 
periodicity angle 9 of the fermion fields. In the degenerate case, i = j, the correlators 
assume the standard form as introduced earlier [24,25], and niij just reduces to the current 
quark mass of a single quark flavour that is prepared by a corresponding choice of equal 
values for the associated hopping parameters, Also the precise definition of the 

lattice derivatives in eq. ()3.3|) matters. As it is written there, do = |(9o + <9o) denotes the 
average of the usual forward and backward derivatives, but as in Ref. [30] we have employed 
the improved derivatives 

do - do (1 - | a 2 d* d ) , d* d - d* Q d (1 - i a 2 d* d ) (3.4) 

as well, which (when acting on smooth functions) have 0(g 2 a 2 ,a 4: ) errors only. 

To enable their numerical calculation, the desired coefficients 6a — bp, b m and the finite 
factor Z = Z m Zp/Zx have to be isolated. This can be achieved by virtue of the identity 



m,i = Z 



~ (m qi j + m q j) + \ b m {am\i + am\j) 

- 3 (6a - bp) a (m q>i + m q j) 2 + O (a 2 ) , (3.5) 



4 This expression for the PCAC masses is only 0(a) improved up to a factor 1 + 5 (6a — 6p)(om qi i + am,j) 



for quark mass dependent cutoff effects. 
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which is obtained if one equates the two available expressions for the 0(a) improved renor- 
malized quark masses, i.e. in terms of the current quark and the subtracted bare quark 
masses (as they appear for the degenerate case in eq. (j2.15[l and thereafter). Forming ra- 
tios of suitable combinations of degenerate and non-degenerate current quark masses in the 
representation (|3.5j) allows to derive the following direct estimators for &a — 6p, o m and Z [30]: 

2 (2mi 2 - m u - m 22 ) , , , n / . s ,„ R \ 

Rap = 7 77 r = o A - bp + O am q i + om qi2 ) , (3.6) 

(mn - m 2 2)(am qi i - om q , 2 ) 

p 4(mi 2 - m 33 ) 

% = t 77 r = Om + O am q ,i + am q;2 , (3.7) 

(mn - m22)(am q ,i - am q ,2) 

with mo,3 = 5(^0,1 + m o,2), apart from other quark mass independent lattice artifacts of 
0(a). For the renormalization constant Z an analogous expression holds even up to 0(a 2 ) 
corrections, 

Rz = — — — — ^- + Q>A - bp - b m )(amu + am 22 ) = Z + O (a 2 ) , (3.8) 
m q ,i - m q , 2 

provided that the correct value for 6a — bp — b m = Rap — Rm (only involving correlation 
functions with mass degenerate quarks) is inserted. 



set 






6 A - 6p 


b m 


Z 


A 


0.132728 


0.4778(7) 


-0.0008(14) 


-0.6217(17) 


1.0941(3) 


B 


0.132711 


0.4621(7) 


-0.0059(22) 


-0.6218(27) 


1.0916(3) 


C 




0.4572(6) 


-0.0057(23) 


-0.6228(28) 


1.0900(3) 


D 


0.132553 


0.4539(5) 


-0.0072(21) 


-0.6260(27) 


1.0882(2) 


E 


0.132395 


0.4508(5) 


-0.0077(25) 


-0.6312(32) 


1.0859(2) 



Table 2: Numerical results on the improvement coefficients 6a — bp and b m and on the renormalization 
constant Z, based on statistics varying between O(900) measurements (A) and O(200) measurements (E). 
These numbers refer to 'choice 1' in eq. 13. 9L where the heavy quark mass is kept at Lmi, ps 0.5, while Lm\ 
indeed turned out to be compatible with zero up to tiny deviations of 0(0.03) (A) and 0(0.01) (sets B - 
E). For C, we interpolated results obtained with L/a = 16, 24 to L/a — 20. Amongst others, 6 m and Z are 
needed to fix the renormalized heavy quark mass in the simulations reported in Section f3. 21 

Concerning the simulation parameters, we argued above that it is most natural to choose 
L/a and (5 exactly as listed in Tabled of the preceding section. Moreover, values for the 
bare quark masses have to be selected. Here we considered two pairs of values for them, 

choice 1 : mo,i = rn c <^ Lm\ ~ , Lm^ i=a 0.5 , 

choice 2 : mo 1 = m c <^ Lm\ ~ , Lm\ 1 ~ 2.6 ; (3-9) 

to comply with the heavy-light notation before we identify the bare quark masses as mo,i = 
7710,1 and 777o, 2 = mo,h with associated PCAC masses mn = mi and m 22 = m^. With the 
first choice (corresponding to z ~ 1) Lm^ is close to the value used in Ref. [30], which 
has the advantage that there this condition was also investigated in perturbation theory 
and the encountered O(a) ambiguities were not enormous. On the contrary, the second one 
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(corresponding to z ~ 5) is motivated by the typical quark mass region that we have to deal 
with when probing our non-degenerate, heavy-light quark system in the next subsection. 
Maintaining these conditions to sufficient accuracy requires some prior simulations to tune 
the hopping parameter belonging to the heavy flavour's bare mass appropriately, whereas 
by K\ = k c taken over from the known values of Ref. [27] the light quark mass is set to zero. 
The final simulation parameters to determine the quantities Rx, X = AP, m, Z, in question 
can be drawn from the triples (L/a, (3, k\) of Tabled specifying the light quark sector (except 
that, only here, T = 3L/2), while the bare heavy quark mass is set through given in the 
second columns of Tables [2] and El There one can also see that the bare current quark mass 
Lnih has been fixed within a few percent to the values dictated by (|3.9|) to keep physics 
constant. 5 



set 






6 A - 6 P 


6m 


Z 


A 


0.126040 


2.7100(6) 


0.0489(3) 


-0.5401(4) 


1.0855(2) 


B 


0.128028 


2.6112(6) 


0.0239(5) 


-0.5621(7) 


1.0867(2) 


C 




2.6709(5) 


0.0151(6) 


-0.5744(9) 


1.0867(2) 


D 


0.129595 


2.5456(5) 


0.0103(5) 


-0.5811(8) 


1.0859(1) 


E 


0.130246 


2.5035(5) 


0.0051(6) 


-0.5927(9) 


1.0845(1) 



Table 3: The same as in Tabled but for 'choice 2' of the heavy quark mass: Lm^ m 2.6. 



b A -b p 
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• choice 1 
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■ choice 2 
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Figure 2: Non-perturbative results for 6a — bp in the considered region of bare couplings for our two choices 
of quark masses (cf. eq. ^3.9^ 1. together with the corresponding results at larger g from [30] (triangles) and 
the prediction from one-loop perturbation theory. The solid line gives the rational fit function that was quoted 
in [30] to well describe the numerical simulation results obtained there. 

The Monte Carlo simulation details and the technical aspects of the analysis to compute 
the observables (|3.6j) - (|3.8|) from the measured fermionic correlation functions are essentially 

° Similar to the situation in Ref. [30] , this is to a good precision equivalent to keeping fixed the corresponding 
renormalized masses Ltur = LZ\m/Zp, because also over our range of considered couplings the entering 
renormalization constant barely varies. 
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• choice 1 

- ■ choice 2 
-0.7 - A NPB595(2001)44 

- -- 1st order PT 
i i i i I i i i i I i i i i I 

0.7 0.75 0.8 0.85 0.9 0.95 1 

St 

Figure 3: The same as in Figure[5]but for the improvement coefficient 6 m . In view of the leading perturbative 
behaviour that one expects to be approached in the limit g% — > 0, the curvature seen in our results hints at a 
more complicated structure of (unknown) higher-order terms. 




the same as in Ref. [30]. In particular, these secondary quantities have been averaged over 
the central timeslices L/(2a), . . . ,(T — L/2)/a to increase statistics, and their statistical 
errors were estimated by the jackknife method. Another issue that deserves to be mentioned 
is the occurrence of the third mass parameter 7710,3 in the determination of 6 m , which has 
to be set properly to ensure the required (but subtle) cancellation leading to the expression 
(|3.7|) for the estimator R m . As in terms of the hopping parameters this amounts to form the 
combination K3 = 2kiK2/(ki + ^2), the condition ^(mo,i +™o,2) _m o,3 = can be satisfied in 
practice only up to roundoff errors, which might cause a systematic uncertainty in simulations 
with single precision arithmetics. But in contrast to Ref. [30], where this effect had to be 
taken into account carefully, the exact numerical value of K3 was now calculated in double 
precision directly in the simulation program itself so that a rounding error contribution 
(reflecting some possible remnant imperfection in the cancellation via 7710,3) to the error on 
R m can be neglected here. 

Our non-perturbative results on b^—bp, b m and Z, which we obtained from the numerical 
simulation data along these lines using improved derivatives throughout, are also collected 
in Tables [21 and |31 6 As a consequence of the fact that all renormalized quantities have been 
fixed in units of L while the ratio L/L max and thereby also the renormalized SF coupling are 
kept constant, the estimates Rx become smooth functions of the bare coupling, g$ = 6/(3. 
This is well fulfilled in Figures El - EJ where our results are shown in conjunction with 
those at larger couplings of the earlier work [30] together with the one-loop perturbative 
predictions [30,37]. One observes that the estimates for 6a — bp and Z corresponding to 
'choice 1' of quark masses are roughly consistent with the fit functions quoted in Ref. [30] 
for larger couplings, if one assumes them to be even valid far beyond the /3-region where the 
underlying data were actually taken. Nonetheless, this was not obvious from the beginning, 

6 We note in passing that memory and processor-topology restrictions of the APEmille parallel computers 
in use prevented us from directly simulating L/a = 20 lattices at an earlier stage of our work. Therefore, we 
employed 16 3 x 30 and 24 3 x 30 lattices and a linear interpolation in 1/L to arrive at the results for set C in 
Tables Hand El 
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Figure 4: The same as in Figure[5]but for the renormalization constant Z. 



AZ F 

0.008 
0.006 
0.004 
0.002 





Ab m 
-0.02 

-0.04 

-0.06 

-0.08 



o -0.1 

0.002 0.004 0.006 0.02 0.04 0.06 0.08 

(aA)' a/L 

Figure 5: Left: Difference of the results on the renormalization constant Z, obtained from the two heavy 
quark mass choices, versus (a/L) 2 . Right: The same for the improvement coefficient bm where, however, the 
ambiguity inherent in any improvement condition imposed is of O(a). 



since despite the quark mass values being comparable in that case the conditions defining 
the associated lines of constant physics are different. 

For b m the situation is more intricate. Here a naive prolongation of the curve from [30] to 
weak couplings (dottet line in Figure |3J) does not give the right g^-behaviour for this region. 
We infer this from the fact that the difference of the dotted line to our results does not seem 
to be compatible with being of O(a). In other words, if the improvement condition of [30] 
were used also in the region of weaker couplings, the points would look quite differently, 
and we have to conclude that a simple approach to one-loop perturbation theory is not an 
adequate representation for the continuation of the triangles in Figure This underlines 
the importance of using particularly adapted improvement conditions, which may be used 
in the coupling range actually relevant for the desired application. The example of b m thus 
illustrates that our redetermination of the 6-coefficients and the Z-factor indeed eliminates 
a source of uncontrollable error. 
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Figure 6: LoFps(Lo, M, go), where the results on b m , Z from both renormalization conditions corresponding 
to eq. ().'-{. 9fl were used to keep the dimensionless RGI heavy quark mass fixed to z — LqM = 9 while varying 
ft — 6/ go and thus a/L. Assuming quadratic scaling violations and discarding the coarsest lattice in the 
(unconstrained) extrapolations, the continuum limits coincide. A more thorough discussion of the continuum 
extrapolations of the observables Fx, X = PS, av, will follow in Section f3. 21 

On the other hand, the results from 'choice 2' (i.e. with the quark mass being 
larger) fall significantly apart throughout, because the a-effects are generically larger in 
that case. But this only reveals an expected, inevitable property of the procedure applied: 
any other estimate Rx (e.g. stemming from a different choice of renormalization condition) 
may yield a different functional dependence upon Oq, but its differences are again smooth 
functions that must vanish in the continuum limit with a rate proportional to a/L (for 
improvement coefficients) or even (a/L) 2 (for renormalization constants). These intrinsic 
0(a n ) ambiguities, n = 1, 2, imply that rather than a numerical value at some given (5, the 
essential information lies in the correct g$ -dependence of the results for the estimators Rx, 
X = AP, m, Z, resulting from working at fixed physics while varying (5. 

To demonstrate the last statement, we also investigated a few alternative improvement 
conditions, which are either realized by defining the estimators Rx with standard instead 
of improved derivatives (as in [30]) or directly by the two different quark mass settings that 
we already have at our disposal through the choices in eq. ().'-{, 9[) . As an example we plot in 
the left part of Figure the difference AZ(gg) = Z(g$)\ choicel - Z(g$) | choice 2 versus (a/L) 2 , 
which clearly shows a linear approach towards zero. Other cases behave similarly, e.g. the 
O(a) ambiguities for Ab m (g^) = b m (a 2 ) | ch oice l - MfluOl choice2 in tlie ri g ht P art of Figure El 
are found to be very small and rapidly decreasing in magnitude as a/L — > 0. 

As a further (and less direct) check for the universality of the continuum limit we con- 
sider a physical quantity that depends on 6 m and Z in a more implict way, namely the 
energy LqTps(Lq, M, go) introduced in eq. (j2.fi |) of the previous section. In fixing z while 
computing -LoFps f° r various lattice resolutions, the just determined results on 6 m and Z 
enter via eq. (|2.18jl . and hence it is interesting to confront the lattice spacing dependences 
of L Q r P s | b m ,z-. choice l and LoFpsU m ,Z:choicc2 with eachother. We did this exercise for z = 9 
(where, due to the large quark mass, a-effects are already very pronounced) and anticipate 
results from the following section in Figure El to display the two data sets and its continuum 
extrapolations linear in (a/L) 2 . From the nice agreement of the continuum limits 7 we infer 

7 The fact that in this case the cutoff effects in I/oTps are larger with b m ,Z from 'choice F (where the 
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once more that our results (on b m and Z) correctly model — within each choice of improve- 
ment condition separately — the respective ^"dependences, entailing convergence to the 
continuum limit with leading corrections of 0(a 2 ). 

3.2 Heavy-light meson energies 

Supposing the parameters (3 and k\ at each L/a (= Lo/a) to be appropriately fixed to comply 
with the conditions g 2 (Lo/2) = 1.8811 = constant and m\ = 0, eqs. 1)2. 13J) and ()2.10j) . we 
still have to prescribe a sequence of dimensionless quark mass values z in order to be able 
to map out the heavy quark mass dependence of the observables Tx, X = PS, av, over a 
reasonable range that encloses the RGI mass scale of the b-quark itself. To this end we split, 
according to the discussion around eq. 1)2.18)) in Section f2. 31 z = LqM into the product 

Z Z A 

z = Lq x h(L ) x Z m x m qh (l + b m am qh ) , Z m = — — . (3.10) 

Z P 

For 6 m (<7g) and Z{g§) we decided to use the results of Table |2] from 'choice 1' of quark 
mass settings, mainly because especially for b m the g$ -dependence is weaker and with the 
corresponding value of the (heavy) quark mass fixed in their determination we more resemble 
the condition that in Ref. [30] was found favourable also from the perturbative point of view. 
Za(^o) i s known in 0(a) improved quenched QCD via the formula [38] 

1 - 0.8496 c/ 2 +0.0610^ „ 

Zx = — o — — , gl < 1 , (3.11) 

1 - 0.7332 c/ 2 ' y ° ~ v ; 

while the required values of Zp (go , Lo/a) have already been quoted in the last column of 
Tableland, as was said after eq. 1)2.20)1 . the universal factor 

h(L ) = 1.531(14) , L = L max /2 , (3.12) 

could be extracted from the data published in Ref. [27]. Given some value of z one is aiming 
at, the relation between the subtracted bare (heavy) quark mass and the hopping parameter, 

am q ,h = \{- -J , (3-13) 

2 V% k c J 

then allows to straightforwardly solve eq. 1)3.10)1 for and, together with K\ as quoted 
in Table ^ of Section 12.31 yields the pairs (/q, Kh) °f hopping parameters, for which the 
numerical simulations with the computation of the heavy-light correlation functions have 
been performed. 8 

Of course, the uncertainties to be associated with the various pieces entering eq. 1)3.10)1 
translate into an error on z, which has to be taken into account also for any quantity regarded 
as a function of z. More precisely, the resulting error on z consists of a go~dependent part 
and a universal, <7o~fnclependent one: while the former comes from the uncertainties of 6 m , Z 

quark mass fixed is smaller) is not so surprising, since via 'choice 2' as improvement condition (where z ~ 5) 
one is closer to the line in parameter space with z = 9 along which LoTps is computed. 

8 Due to the fact that the non-perturbative values for & m in the relevant /3-range lie around —0.6, the 
relation z = Loh(Lo)Z ul m qt h{l + 6 m cwTiq,h) can not be inverted in favour of am q ,h (and thereby Kh) for 
arbitrarily high z-values. In case of our largest z (= 13.5), for instance, this already restricts the possible 
inverse lattice spacings in Table |l]to L/a = 16 — 32. 
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and Zp quoted in the tables plus an error (of 0.8% at f3 « 7.4 down to 0.4% for (3 > 7.8 [38]) 
on Za, the latter stems from the overall uncertainty of 0.9% on h(Lo) in the continuum 
limit and hence has only to be added in quadrature after a continuum extrapolation of the 
respective z-dependent quantity under study. 9 

The final parameter sets employed in the Monte Carlos simulations consist of the triples 
(L/a, (3, k\) in Tableland the Kh-values we arrived at as described before; they are listed 
together with the corresponding values of the dimensionless RGI heavy quark mass, z, and 
the ^"dependent error part of the latter in Table |1] at the end of this section. For the 
technical details of the runs to produce the numerical data on the SF heavy-light meson cor- 
relators, from which the logarithmic derivatives (|2.6|) and (|2.7fl are evaluated, we refer to the 
simulations reported in Ref. [31] (and its Appendix A. 2. in particular) to non-perturbatively 
renormalize the static-light axial current. The number of measurements in the statistical 
samples is comparable to what was accumulated to get the results in Tables [2] and El 

Table |1] also contains the numerical results for LoFps and LoT av at all values of z that 
were considered. To ease notation, we set 

n x (u,z,a/L)=L T x (L ,M,g ) , X = PS,av, (3.14) 

g 2 (L /2)=u,L M=z 

and denote its continuum limits as 

ujx(u, z) = lim Qx(u,z,a/L), X = PS,av. (3.15) 

a/L->0 

In Table E] both the statistical errors of fix (in square brackets) as well as the combined 
statistical and go~dependent, z-induced uncertainties (in parentheses) are given, the latter 
being obtained by including the propagation of the ^"dependent part of the error on z. 
As an outcome of this error analysis we find that summing up the various errors on the 
factors in eq. (|3.10j) quadratically (since they are uncorrelated) and multiplying with the 
numerically estimated slopes | [ 9f2x/<9.z ] (9o)\ yields a contribution of about 0.3% or less to 
the final uncertainties on fix- 
As we work in the 0(a) improved theory, the numbers for Q-ps and fl av can now be 
extrapolated linearly in a 2 to the continuum limit. 10 However, based on the experience made 
in perturbation theory (see Section 5 of Ref. [39]) that the discretization errors primarily 
depend on the mass of the heavy quark, aM = z x a/L, we expect some deviation from a 
leading linear behaviour in a 2 for the coarsest lattices as z grows. Since our non-perturbative 
data qualitatively confirm this picture, we again adopt perturbation theory as a guide to get a 
rough estimate for the quark mass value where this deviation sets in. In [39] the heavy quark 
mass dependence of the discretization errors in a typical matrix element of the heavy-light 
axial current, similarly constructed from SF correlators as our energy observables, indicates 
a breakdown of O(a) improvement beyond (am^ s ) 2 ~ 0.2, which with m^ s ~ 0.7M [40] 
approximately corresponds to aM ~ 0.64. From this we deduce the following two-step 
criterium to carefully perform the continuum extrapolations of fix, X = PS,av, for the 
various values of z: 



9 Note that also this error can in principle be reduced further by increasing the precision of the continuum 
step scaling functions of Ref. [27] . 

10 As in [27,31], the influence of the only perturbatively known SF-specific boundary improvement coeffi- 
cients c t and c t is negligible at the level of our precision. 
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1. In view of the ^-dependence of the size of the lattice artifacts, only allow for fits which 
start at and beyond the minimal L/a such that aM < 0.6 approximately holds. 

2. Among these, the n-point fit must agree with the (n + l)-point one within errors but 
the former, which omits the coarsest of the lattices meeting 1. and thereby has the 
larger error, gives the final estimate of the continuum limit. 




Figure 7: Continuum extrapolations of f2ps(l. 8811, z, a/L) (left) and fi av (1.8811, z, a/L) (right) for rep- 
resentative values of z. As explained in the text, the number of coarsest lattices to be skipped in the fits 
depends on z. The final continuum limits are taken from the solid fit functions, which are compatible with 
the dashed extrapolations omitting one point less. 

The resulting continuum extrapolations for the subset z £ {5.15,6.0,9.0,13.5} of avail- 
able z -values are displayed in Figure [Tj As it turns out when applying this fitting procedure, 
the continuum limits then have to be extracted from fits discarding the L/a = 12 lattice in 
case of z € {3.0,3.8,5.15,6.0,6.6}, the L/a = 12,16 lattices in case of z = 9.0 and even the 
L/a = 12, 16,20 lattices for the heaviest mass, z = 13.5. The numbers of Table0]in Italics 
are the continuum limit results for all z, based on these fits. 

To further corroborate our results from this prescription, we in addition considered an 
alternative fit ansatz that also accounts for a term cubic in a. Again, it is to some extent 
guided by the aforementioned finding of Ref. [39] that in perturbation theory the cutoff 
effects in the regime of large quark masses are approximately a function of aM and not of 
aM and a/L separately. In the O(a) improved case at hand (and here for the example of 
O av , omitting the argument u for a moment), such an ansatz may therefore be written as: 



f2 av (z,a/L) = uj av (z) 



where by a z-dependent parameter pi we after all admit a still more general form for the 
a 2 -term. Performing a simultaneous fit of all available f2 av data (no cut on aM), we then 
arrive at continuum limits in complete accordance with the former, having comparable or 
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even smaller errors: uj^(z) e {2.86(1), 3.34(1), 4.15(3), 4.65(2), 5.00(2), 6.39(3), 8.99(6)} for 
z = 3.0, . . . , 13.5. But keeping in mind that for aM ~ 0.6 — 0.8 the perturbative a-expansion 
entirely breaks down [39], we take the continuum limits of Table 0] from the linear extrapola- 
tions in (a/L) 2 imposing a cut on aM as a safeguard against any uncontrollable higher-order 
behaviour as our final results, which within their (larger) errors are moreover also consistent 
with the values at smallest lattice resolution. All told, we thus are confident that possible 
systematic uncertainties in the fitting procedure are already well covered by these estimates. 

For future reference (i.e. in particular for the non-perturbative determination of the fa- 
quark's mass in Ref. [9]) we introduce for the dimensionless, spin-averaged heavy-light meson 
energy in small volume the further abbreviation 

f2(it, z, a/L) = O av (ii, z, a/L) , u>(u, z) = w av (u, z) . (3-17) 

To parametrize the z-dependence of our continuum values w(1.8811, z) by a smooth fit 
function, we make the ansatz a^z + ai + 02/2, which is justified by the theoretical expectation 
that heavy-light meson correlation functions of the type studied here decay with a leading 
term proportional to the heavy quark mass, up to some low-energy scale of O(Aqcd) and 
1/m-corrections. Since the simulation data for the various z at given (3 were produced on 
the same gauge field backgrounds, possible correlations in z had to be taken into account in 
this fit. Hence we performed it on basis of the jackknife samples that were built from the 
raw data and passed through the whole analysis, and we end up with the parametrization 

uj(lMll,z) =a z + a 1 + a 2 ~, a = 0.581 , ai = 1.226 , a 2 = -0.358, (3.18) 

z 

the graph of which is shown in the upper diagram of Figure[S]to well represent our data. In the 
10 I 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 _| 




2 4 6 8 10 12 14 

z 

Figure 8: Continuum limit values w(1.8811,,z) of the spin-averaged energies £1(1. 8811, z, a/L), at fixed 
coupling g 2 (L /2) = 1.8811, as a function of z — LoM and its fit function (top). The bottom graph has the 
linear term, a^z, subtracted and also includes the pseudoscalar case, wps(1.8811, z), for comparison. 

interval 5.2 < z < 6.6, which for instance is the relevant z-range to calculate the RGI b-quark 
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mass by means of the non-perturbative matching to HQET in Ref. [9], this parametrization 
describes w(1.8811, z) with a precision of about 0.5%. As already mentioned, a further 
uncertainty of 0.9% that enters only indirectly through the argument z of the function to 
(and thus remains to be added in quadrature at the end) originates from the universal part 
h(Lo) in eq. (j3.10|) of the renormalization factor relating the bare to the RGI quark mass. 

For later use of the result l|3.18|) it is also necessary to have a numerical estimate for 
the derivative u/(1.8811, z) w.r.t. u, in order to compensate for a possible slight mismatch 
in the imposed condition of fixed renormalized coupling. From an additional simulation 
with L/a = 24 and a nearby coupling of <? 2 (Lo/2) = 1.95 (and assuming the a-effect of the 
derivative to be negligible) we found this derivative to be constant in the central region 



6.0<,<6.6: 

ou 



0.70(1) 



(3.19) 



m=1.8811 



which is the required one in the context of the immediate application discussed in [9]. 

Finally, we give as well the fit result for ups (1.8811, z) for comparison, i.e. if in the 
■ 14|) the effective energy in the pseudoscalar channel alone, Tps(Lq, M, go), would 



definition 

be used instead of the spin-averaged combination. By an analogous analysis with the same 
ansatz for the fit function as in eq. (|3.18|) we obtain the coefficients 



a = 0.587 , ai = 1.121 , a 2 = -1.306 (for T x = r PS ) . 



(3.20) 



The fact that the absolute value of a 2 is evidently larger than in the spin-averaged case is 
illustrated by the shape of the subtracted fits in the bottom diagram of Figure[S]and confirms 
our earlier claim that with Tx 



T av the matching condition (|2.11j) becomes independent of 
the coefficient of the chromomagnetic field strength term in the HQET Lagrangian also at 
first order inl/m. Therefore, the use of Q = O av = LoT av (and eq. (|3.18|) for its heavy quark 
mass dependence in the continuum limit) is preferable for later applications to HQET, since 
the 1/m-errors arising at leading order will generically be smaller then. 
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z 


set 


K h 




f2 PS 




3.0 


A 


0.129571 


3.000(12) 


2.534(11)[5] 


2.916(11)[6] 




B 


0.130445 


3.000(13) 


2.485(10) [4] 


2.869(11)[6] 




C 


0.130876 


3.000(13) 


2.465(10) [5] 


2.856(12)[8] 




D 


0.131115 


3.000(13) 


2.471(11)[7] 


2.863(14) [12] 




E 


0.131349 


3.000(13) 


2.458(13) [10] 


2.858(18)[16] 




CL 






2451(13) 


2.852(18) 


3.8 


A 


0.128310 


3.800(16) 


3.094(13) [5] 


3.408(14) [6] 




B 


0.129552 


3.800(16) 


3.039(12) [4] 


3.356(13)[6] 




C 


0.130185 


3.800(16) 


3.019(12) [5] 


3.341(13)[8] 




D 


0.130553 


3.800(16) 


3.022(12) [7] 


3.345(15)[11] 




E 


0.130940 


3.800(17) 


3.011(15)[10] 


3.342(19)[16] 




CL 






3.002(16) 


3.334(19) 


5.15 


A 


0.126055 


5.150(21) 


4.025(17)[5] 


4.262(17)[6] 




B 


0.127991 


5.150(22) 


3.948(16)[5] 


4. 188(16) [6] 




C 


0.128989 


5.150(22) 


3.920(15) [5] 


4.166(16) [8] 




D 


0.129586 


5.150(22) 


3.917(16) [7] 


4.163(18)[11] 




E 


0.130242 


5.150(23) 


3.905(18)[11] 


4.161(21)[16] 




CL 






3.890(20) 


4.146(22) 


6.0 


A 


0.124528 


6.000(25) 


4.619(20)[5] 


4.821(20)[6] 




B 


0.126967 


6.000(25) 


4.513(18)[5] 


4.720(18) [6] 




C 


0.128214 


6.000(25) 


4.476(17) [5] 


4.688(18) [8] 




D 


0.128964 


6.000(25) 


4.467(18) [8] 


4.679(19)[11] 




E 


0.129796 


6.000(27) 


4.453(20)[11] 


4.674(23) [15] 




CL 






4430(22) 


4.652(24) 


6.6 


A 


0.123383 


6.600(27) 


5.049(22) [5] 


5.232(22)[6] 




B 


0.126222 


6.600(28) 


4.912(20) [5] 


5. 100(20) [6] 




C 


0.127656 


6.600(28) 


4.866(19)[5] 


5.059(20) [8] 




D 


0.128518 


6.600(28) 


4.852(19)[8] 


5.045(21)[10] 




E 


0.129477 


6.600(30) 


4.835(22)[11] 


5.037(24) [15] 




CL 






4.806(24) 


5.009(26) 


9.0 


A 


0.117762 


9.000(39) 


7.022(31)[5] 


7.146(31)[5] 




B 


0.122987 


9.000(39) 


6.554(27) [5] 


6.690(27) [6] 




C 


0.125309 


9.000(38) 


6.429(25) [6] 


6.572(26) [7] 




D 


0.126670 


9.000(38) 


6.378(26) [8] 


6.521(26)[10] 




E 


0.128175 


9.000(40) 


6.337(28) [12] 


6.487(29) [15] 




CL 






6.277(47) 


6.429(48) 


13.5 


B 


0.113764 


13.500(69) 


10.721(47)[5] 


10.801(47)[5] 




C 


0.120094 


13.500(60) 


9.603(40) [6] 


9.695(40) [7] 




D 


0.122808 


13.500(58) 


9.324(38) [8] 


9.418(38) [9] 




E 


0.125575 


13.500(61) 


9.136(40)[12] 


9.235(41)[14] 




CL 






8.89(10) 


9.00(11) 



Table 4: Heavy quark masses z = LoM with go-dependent error and associated results for flx(u,z,a/L), 
X = PS, av, u — 1.8811, with the total go-dependent part of the error in parentheses and the only statistical 
one of fix in square brackets. Continuum limits (see text) are displayed in Italics. 
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4 Conclusions 



In this work we have computed the non-perturbative heavy quark mass dependence of ef- 
fective energies derived from heavy-light meson correlation functions by means of numerical 
simulations. A particular aspect of this computation is the use of a physically small volume 
(of a linear extent of 0(0.2 fin)): it is a prerequisite to treat the heavy quark flavour in lattice 
regularization as a relativistic particle and thus may be looked at as a 'device' that serves 
to non-perturbatively renormalize HQET with the method of matching the effective theory 
to small- volume QCD, proposed and applied in Ref. [9]. 

Among the several improvement coefficients and renormalization constants that are 
needed to accurately keep fixed the renormalized (heavy) quark mass while the observables 
of interest approach their continuum limits, we determined those relating the renormalized 
to the subtracted bare quark mass in a range of weaker couplings relevant here (7.4 < 
(3 = 6/c/q < 8.2) with high precision. Even in the region of the b-quark mass and slightly 
beyond we still find our results on the meson energies under study to be rather mildly cutoff 
dependent so that the continuum extrapolations can be well controlled. 

The quenched numerical results obtained in Section 13.21 directly pass into the — for the 
first time entirely non-perturbative — calculation of the b-quark mass in the static approxi- 
mation of [9,41]. Moreover, our quantitative investigation of the quark mass dependence in 
relativistic heavy-light systems in finite volume can easily be extended to other heavy-light 
bilinears and matrix elements, which then opens the possibility to perform genuinely non- 
perturbative tests of HQET and to estimate the size of the 1/m-corrections to the static 
theory. We will focus on these issues in a separate publication [12]. 

Finally, considering the small-volume investigation presented here in the light of the 
more general framework [9] of a non-perturbative matching between HQET and QCD, it is 
an obvious practical advantage that they can be transfered to also include dynamical quarks 
in the (hopefully near) future without requiring exceedingly large computing resources. 
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